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We study the classical problem of planar shock refraction at an oblique density disconti- 
nuity, separating two gases at rest. When the shock impinges on the density discontinuity, 
it refracts and in the hydrodynamical case 3 signals arise. Regular refraction means that 
these signals meet at a single point, called the triple point. 

After reflection from the top wall, the contact discontinuity becomes unstable due to 
local Kelvin-Helmholtz instability, causing the contact surface to roll up and develop 
the Richtmyer-Meshkov instability. We present an exact Riemann solver based solution 
strategy to describe the initial self similar refraction phase, by which we can quantify the 
vorticity deposited on the contact interface. We investigate the effect of a perpendicu- 
lar magnetic field and quantify how addition of a perpendicular magnetic field increases 
the deposition of vorticity on the contact interface slightly under constant Atwood Num- 
ber. We predict wave pattern transitions, in agreement with experiments, von Neumann 
shock refraction theory, and numerical simulations performed with the grid-adaptive code 
AMRVAC. These simulations also describe the later phase of the Richtmyer-Meshkov in- 
stability. 



1. Introduction 

We study the classica l problem of regular r efraction of a shock at an oblique density 
discontinuity. Long ago, von Neumann (|l943h deduced the critical angles for regularity 
of the ref raction, wh ile Taubl ( 1947 ) found relations between the angles of refraction. 
Later on, [Henderson (|l966l) extended this work to irregular refraction by use of pola r 
diagrams. An example of an early shock tube experiment wa s performed bv lJahn ( 19561 ). 
Amongst many others, Abd-El-Fattah fc Henderson (197830) performed experiments in 
which also irregular refraction occured. 

In 1960, Richtmyer performed the linear stability analysis of the interaction of shock 
waves with density discontinuities, and concluded that the shock-ac celerated contact is 
unstable to perturbations of all wavelenghts, for fast-slow interfaces ( Richtmyeil ( 196Cll )). 
In hydrodynamics (HD) an interface is said to be fast-slow if 77 > 1, and slow-fast oth- 
erwise, where r] is the density ratio across the interface (figure [T]). The instability is not 
a classical fluid instability in the sense that the perturbations g row linear l y and not 
exponentially. The first experimental validation was performed by iMeshkovl ( 1969t ). On 
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the other hand, according to linear analysis the interface remains stable for slow-fast 
interfaces. This misleading result is only valid in the lin ear phase of the process and 
near th e triple point: a wid e range of experimenta l (e.g. lAbd-El-Fattah fc Henderson 
( 197861) ) and numerical (e.g. iNouragliev et all results show that also in this case 

the interface becomes unstable. The growth rates obtained by linea r theory compare 
poorly to experimentally determined growth rates (jSturtevanti The governing 

instability is referred to as the Richtmyer-Meshkov inst ability (RMI) an d is nowadays a 
topi c of research in e . g. ine rtial confinement fusion ( e.g. lOron et all (|l999i) ). astrophysics 



e.g. Kifonidis et al. ( 20061)), an d it is a common test problem for numerical codes ( e.g. 



van der Hoist fc Keppe^(|2007i )). 



In essence, the RMI is a local Kelvin-Helmholtz in s tabilit y, due to the deposition 
of vorticity on the shocked contact. Hawlev fc Zabuskv ( 19891 ) formulate an interesting 
vortex paradigm, which desc ribes the p rocess of shock refraction, using vorticity as a 
central concept. Later on. iSamtanev et al.. (1998 ) performed an extensive analysis of the 
baroclinic circulation generation on shocked slow-fast interfaces. 

A wide range of fields where the RMI occurs, involves ionized, quasi-neutral plas- 
mas, where the magnetic field plays an important role. Therefore, more rec ently there 
has be en some research done on the RMI in magnetohydrodynamics (MHD). ISamtanevI 
(|2003 ^ proved by numerical simulations, exploiting Adaptive Mesh Refinement (AMR), 
that the R MI is suppressed in pla nar MHD, when the initial magnetic field is normal to 
the shock. IWheatlev et al. I (l2005l) solved the problem of planar shock refraction analyt- 
ically, making initial guesses for the refracted angles. The basic idea is that ideal MHD 
does not allow for a jump in tangential velocity, if the ma gnetic field component norm al 
to the contact discontinuity (CD), does not vanish (see e.g. Goedbloed fc Poedtd ( 2004I )). 
The solutio n of the Riemann problem in ideal MHD is well-studied in the literature (e.g. 
Lax ( 1957[ )). and due to the existence of three (slow, Alfven, fast) wave signals instead 



of one (sound) signal, it is much richer than the HD case. The Riemann problem usually 
considers the self similar temporal evolution of an initial discontinuity, while we will con- 
sider stationary two dimensional conditions. The interaction of small pertu rbations with 
MHD (switch-on a nd switch-off) shocks w as studied both analytically by iToddl (|l965h 
and numerically by Chu fc Taussig! (1967). Later on, the evolutionarity of intermediate 
shocks, which cross the Alfven speed, has been studied extensively. I ntermediate shocks 
are unstable under small p erturbations, and are thus not evolutionary. iBrio fc Wu ( 1988h 
and lDe Sterck et ali (|l998l ) found intermediate shocks in respectively one and two dimen- 
sional simulations. Th e evolutionary condition became controversial and amongst others 
Mvong fc Roe (I1997q||6I ) argu e that the evolutionary condition is not relevant in dissipa- 
tive MHD. IChao et aLl ( iggs') reported a 2 ^ 4 intermediate shock observed by Voyager 
1 in 1980 and iFeng fc Wan g (2008) recognised a 2 - -» 3 intermediate sho ck, which was 
observed by Voyager 2 in 1979. On the other hand, Barmin et al. ( 1996h argue that if 
the full set of MHD equations is used to solve planar MHD, a small tangential distur- 
bance on the magnetic field vector splits the rotational jump from the compound wave, 
transforming it into a slow shock. They investigate the recons truction process of the non- 
evolu tionary compound wave into evolutionary shocks. Also iFalle & Komissarovl ()1997L 



20011 ) do not reject the evolutionary condition, and develop a shock capturing scheme for 



evolutionary solutions in MHD, However, since all the signals in this paper are essentially 
hydrodynamical, we do not have to worry about evolutionarity for the setup considered 
here. 

In this paper, we solve the problem of regular shock refraction exactly, by developing a 
stationary two-dimensional Riemann solver. Since a normal component of the magnetic 
field suppresses the RMI, we investigate the effect of a perpendicular magnetic field. 
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Figure 1. Initial configuration: a shock moves with shock speed M to an inclined density 
discontinuity. Both the upper and lower boundary are solid walls, while the left and the right 
boundaries are open. 



The transition from slow-fast to fast-slow refraction is described in a natural way and 
the method can predict wave pattern tr ansitions. We also perform nume ri cal simulations 



using the grid- adaptive code AMRVAC ( van der Hoist fc Keppend ( 2007 ): Keppens et al 
((200ai l. 

In section 2, we formulate the problem and introduce the governing MHD equations. 
In section 3, we present our Riemann solver based solution strategy and in section 4, 
more details on the numerical implementation are described. Finally, in section 5, we 
present our results, including a case study, the prediction of wave pattern transitions, 
comparison to experiments and numerical simulations, and the effect of a perpendicular 
magnetic field on the stability of the CD. 



2. Configuration and governing equations 

2.1. Problem setup 

As indicated in figure [l] the hydrodynamical problem of regular shock refraction is 
parametrised by 5 independent initial parameters: the angle a between the shock normal 
and the initial density discontinuity CD, the sonic Mach number M of the impinging 
shock, the density ratio 77 across the CD and the ratios of specific heat 7; and 7^ on both 
sides of the CD. The shock refracts in 3 signals: a reflected signal (R), a transmitted 
signal (T) and a shocked contact discontinuity (CD), where we allow both R and T to 
be expansion fans or shocks. Adding a perpendicular magnetic field, B, also introduces 
the plasma-/? in the pre-shock region, 

= (2-1) 

which is in our setup a sixth independent parameter. As argued later, the shock then still 
refracts in 3 signals (see figure [l]): a reflected signal (R), a transmitted signal (T) and a 
shocked contact discontinuity (CD), where we allow both R and T to be expansion fans 
or shocks. 



2.2. Stationary MHD equations 

In order to describe the dynamical behaviour of ionized, quasi-neutral plasmas, we use 
the framework of ideal MHD. We thereby neglect viscosity and resistivity, and suppose 
that the length scales of interest are much larger than the Debye leng th and there are 
enough particles in a Debye sphere (see e.g. lCoedbloed fc Poedtd ()2004l )). As written out 
in conservative form and for our planar problem, the stationary MHD equations are 
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Figure 2. Left: A stationary shock, seperating two constant states across an inclined planar 
discontinuity. Right: The eigenvalues of the matrix A from (|3.12|l correspond to the refracted 
signals. 



dx dy 
where we introduced the flux terms 



d d 

F + — G = 0, (2.2) 



F = 

and 

G 



t 



pVx,pvl+p+^,pv^Vy,v^{—!-jP + p ^ ^ + B^),v^B,v^-fp\ , (2.3) 



t 



pVy,pV^Vy,pV^+P+ —,Vy{-—^P + P ^ ^ + ) , Vy B , VyJ P ^ . (2.4) 



The applied magnetic field B — (0, 0, B) is assumed purely perpendicular to the flow and 
the velocity v = {vx,Vy,0). Note that the ratio of specific heats, 7, is interpreted as a 
variable, rather than as an equation parameter, which is done to treat gases and plasmas 
in a simple analytical and numerical way. The latter equation of the system expresses 
that V • (7pv) = 0. Also note that V • B = is trivially satisfied. 

2.3. Planar stationary Rankine-Hugoniot condition 

We allow weak solutions of the system, which are solutions of the integral form of the 
MHD equations. The shock occuring in the problem setup, as well as those that later on 
may appear as i? or T signals obey the Rankine-Hugoniot conditions. In the case of two 
dimensional stationary flows (see figure H]), where the shock speed s = 0, the Rankine- 
Hugoniot conditions follow from equation (|2.2|) . When considering a thin continuous 
transition layer in between the two regions, with thickness S, solutions of the integral 
form of equation (|2.2p should satisfy lim (^F + ■^G)dl — 0. For vanishing thickness 
of the transition layer this yields the Rankine-Hugoniot conditions as 

- lim (^^F - -^^g] dl = (2.5) 
"5^0 ./i V sm 6 dl cos ch dl J 
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Figure 3. The wave pattern during interaction of the shock with the CD. The upper and 
lower boundaries are rigid walls, while the left and right boundaries are open. 



where ^ — tan (p and is the angle between the x-axis and the shock as indicated in 
figure m The symbol [[ ]] indicates the jump across the interface. 

3. Riemann Solver based solution strategy 

3.1. Dimensionless representation 

In this section we present how we initialise the problem in a dimensionless manner. In 
the initial refraction phase, the shock wil introduce 3 wave signals (R, CD, T), and 
2 new constant states develop, as schematically shown in figure [31 We choose a rep- 
resentation in which the initial shock speed s equals its sonic Mach number M. We 
determine the value of the primitive variables in the post-shock region by applying 
the stationary Rankine-Hugoniot conditions in the shock rest frame. In absence of a 
magnetic field, we use a slightly different way to nondimensionalise the problem. Note 
Ui — {pi,Vx,i,Vy^i,ptotA, Bi,"fi), where the index i refers to the value taken in the i—th. 
region (figure [3]) and the total pressure 



In the HD case, we define p ~ I and p = 7; in Ui. Now all velocity components are scaled 
with respect to the sound speed in this region between the impinging shock and the initial 
CD. Since this region is initially at rest, the sonic Mach number Af of the shock equals its 
shock speed s. When the shock intersects the CD, the triple point follows the unshocked 
contact slip line. It does so at a speed Vtp — (M,M tan a), therefore we will solve the 
problem in the frame of the stationary triple point. We will look for selfsimilar solutions 
in this frame, u = u((/)), where all signals are stationary. We now have that = Vx — M 
and Vy = Vy — Af tana, where v refers to this new frame. From now on we will drop the 
tilde and only use this new frame. We now have Ui = (7;, — Af, — Af tana, 1, 0, 7;)* and 
U5 = {rjji, —M, —Mtana, 1,0, jrY . The Rankine-Hugoniot relations now immediately 
give a unique solution for U2, namely 



t 

[[F]] = U[G]] 



(2.6) 



52 

Ptot =P + —■ 



(3.1) 



U2 = 



( 



(7; - + 2 

ill + l)M 



, —Mtana, 



11 + 1 



,0,7/ 



) 



(3.2) 



(7i-l)Af2 + 2' 



11 + 1 
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In MHD, we nondimensionalise by definining B = 1 and p — in region 1. Again 
all velocity components are scaled with respect to the sound speed in this region. We 

now have that Ui ~ (^^i —M, —Mtana, ^y^, 1, 7/^ and from the definition of j], U5 = 
^ — Af, —Mtana, ^-j^, 1, 7r^ ■ The Rankine-Hugoniot relations now give the follow- 



ing non-trivial solutions for U2: 

U2 = ( - 

where 



2u> 



Lu, — Mtana, p2 



2^' 



-M 



-,7/ 



Auj + B 



P2 



(3.3) 
(3.4) 



Clu + D' 

is the thermal pressure in the post shock region. We introduced the coefficients 

A = ^i {(3\A-/fM^ - 2^iM^ -^i-i) + (3 ((^2 ^ _ ^^^,^2 _ 3) - 7; 2) , (3.5) 

B = (7, - 1)M (/3(M2(7f + 77O - 27, + 4) - 27, + 4) , (3.6) 

C = 27,(7, + 1) (/3((7i - 1)M' + 2) + 2) , (3.7) 

I? = 4(7, + 1)(7, -2)M. (3.8) 

The quantity 

7K7;-l)/3M^ + 27i(/3 + l)± 



W 



(3.9) 



27/(7/ + 1)/3M 

is the normal post-shock velocity relative to the shock, with 

W = P^M^i-ff - jf) (m2(7/ _ 1) + 4) + /37/(4M2(4 + 7, - jf) + 87,) + 47^(3.10) 

Note that lo must satisfy — M < lu < to represent a genuine right moving shock. We 
choose the solution where lu ~ since the alternative, lu — w_ is a degenerate solution 
in the sense that the hydrodynamical limit lim = 0, which does not represent a 

rightmoving shock. 

3.2. Relations across a contact discontinuity and an expansion fan 
Rewriting equation (j2.2p in quasilinear form leads to 

u, + (F^-i • Gu) Uy = 0. (3.11) 

In the frame moving with the triple point, we are searching for selfsimilar solutions and we 
can introduce £. = — tan0, so that u — u(^). Assuming that ^ 1-^ u(^) is differentiable, 
manipulating (|3.1ip leads to Au^ — ^u^. So the eigenvalues Xi of A represent tan0, 
where (j) is the angle between the refracted signals and the negative x-axis. The matrix 
A is given by 





Vy_ 


pVy 




Vy 1 








\ 


Vx 




— c 


•Ux V^—C' 








V^Vy 




Vy 1 














— ^ T 




^ T 

Vy_ 

fx 


p v'^—c'^ 

PVj: 









pC^Vy 


PC^Ux 


V^Vy 












vi—c'^ 




f X 










BVy 


SDx 


Vy B 1 


Vy_ 




Vy_ 

fx 




I 








•Ux P 'Ux~C^ 




Vx 




/ 



(3.12) 



Regular shock refraction 



and its eigenvalues are 



Ai, 



2,3,4,5,6 



{- 



y y y y 



VxVy 



-}, 



where the magnetosonic speed c = -^/u^ + and the sound speed Vg 



(3.13) 



^ and the 

p 



Alfven speed Va 



Since A has 3 different eigenvalues, 3 different signals will arise. 



When exists and 7^ 0, i.e. inside of expansion fans, is proportional to a right 
eigenvector of A. Derivation of ^ = with respect to ^ gives (VuAi) • ua = 1 and thus 
we find the proportionality constant, giving 



(3.14) 



VuAi • ri 

While this result assumed continuous functions, we can also mention relations that 
hold even across discontinuities like the CD. Denoting the ratio 



dx—Xj dy 



K, it follows 

(li • rj)K = Kdi_j, where li and are respectively left and right 



that [li • du 

eigenvectors corresponding to A^. Therefore, if i 7^ j. 



(3.15) 



From these general considerations the following relations hold across the contact or shear 



wave where the ratio -^ = 77- 



VydVx 
VydVx 



VxdVy 
VxdVy 



^ — 

pvi. 

c\/ v-^ —c 



p^i 



■dptot 
■dptot 



(3.16) 



Since v ^ c, otherwise all signals would coincide, it follows immediately that the total 
pressure ptot and the direction of the streamlines ^ remain constant across the shocked 
contact discontinuity. 

These relations across the CD allow to solve the full pr oblem using an iterative proce- 
dure. Inspired by the exact Ricmann solver described in iTorol (|l999l ). we first guess the 
total pressure p* across the CD. i? is a shock when p* is larger than the post-shock total 
pressure and T is a shock, only if p* is larger than the pre-shock total pressure. Note 
that the jump in tangential velocity aross the CD is a function of p* and it must vanish. 
A simple Newton-Raphson iteration on this function [[^]](p*), finds the correct p* . We 
explain further in section 3.5 how we find the functional expression and iterate to even- 
tually quantify <^ij, ^r, 4'CD and the full solution \i{x^y,t). From now on p* represents 
the constant total pressure across the CD. 

Similarly, from the general considerations above, equation (|3.15[) gives that along ^ = 
v^vy±cy/v^-c^ ^T^^ following relations connect two states across expansion fans: 



hdptot = 0, 

dp tot = 0, 



dp - -pi<. 

VxdVx + V 

- pdptot + Ptotpdl + Ptotldp = 0, 
-Bdptot + ilP + dB = 0, 



(3.17) 



(7P- 

VydVx - VxdVy ± 



-dp tot = 0. 



These can be written in a form which we exploit to numerically integrate the solution 
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through expansion fans, namely 



-dptot, 
■dp tot, 



(3.18) 



Pi = Pe + Jp^^^^ ^dptot, 

■ ■ '^Ptot.e pV-^C 

Pi ^Pe + Jp^^^^ ^dptot, 

The indices i and e stand respectively for internal and external, the states at both sides 
of the expansion fans. The upper signs hold for reflected expansion fans (i.e. of type R), 
while the lower sign holds for transmitted expansion fans (i.e. of type T). 



3.3. Relations across a shock 

Since the system is nonlinear and allows for large-amplitude shock waves, the analysis 
given thus far is not sufficient. We must include the possibility of one or both of the R 
and T signals to be solutions of the stationary Rankine-Hugoniot conditions (equation 
(|2.6[) ). The solution is given by 



where 



Pi 



7+1 Ptot,, 



7-1 



Pe, 



7 + 1 Vtot,, 
_ ^+(P*-Ptot,<,) 

^•^ pc(i'm,c5T-"y.<!) ' 

I P'-Ptot.c 



7-1 

7+1 ^ Plot,, 



7-1 



7 + 1 Ptot,< 

0i?/T = atan{C+/_), 



(3.19) 



(3.20) 



and 



(7 - l)Ptot.e + (7 + 1)P* 
2pe 



(3.21) 



Again the indices i and e stand respectively for internal and external, the states at both 
sides of the shocks. The upper signs holds for reflected shocks, while the lower sign holds 
for transmitted shocks. 



3.4. Shock refraction as a Riemann problem 

We are now ready to formulate our iterative solution strategy. Since there exist 2 in- 
variants across the CD, it follows that we can do an iteration, if we are able to express 
one invariant in function of the other. As mentioned earlier, we choose to iterate on 
p* — ptot.3 = PtotA- This is the only state variab;e in the solution, and it controls both 
R and T. We will write (jiR = (j)}i{u2,p*) and (f>T = 0t(u5,p*), U3 = U3(u2,p*) and 
U4 = U4(u5,p*). The other invariant should match too, i.e. — = 0. Since U2 
and U5 only depend on the input parameters, this last expression is a function oi p* . 
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Iteration on p* gives p* and (j)R = (pB.{p*)i 4't — (j)T{p*), U3 = U3(p*) and U4 = U4{p*) 
give (bcD — atan^^^ = atan^^^, which solves the problem. 



3.5. Solution inside of an expansion fan 

The only ingredient not yet fully specified by our description above is how to determine 
the variation through possible expansion fans. This can be done once the solution for p* is 
iteratively found, by integrating equations (j3.18|) till the appropriate value oiptot- Notice 

that the location of the tail of the expansion fan is found by tan{4)taii) — "'^■'"'"''^^'^^'V^"' 



and the position of (j>head is uniquely determined by tan{(f)head) = ^"'^^^^"^^'2^^° 

Inside an expansion fan we know u(ptot), so now we need to find ptot{4')i in order to 
find a solution for u(0). We decompose vectors locally in the normal and tangential 
directions, which are respectively referred to with the indices n and t. We denote taking 
derivatives with respect to as '. Inside of the expansion fans we have some invariants 
given by equations (j3.17p . The fourth of these immediately leads to as an invariant. 
Eliminating ptot from dp— -^dptot — and — Bdptot + {lP+ B^)dB ~ yields the invariant 
■g, and combining these 2 invariants tells us that the entropy S* = ^ is invariant. The 
stationary MHD equations (|2.2p can then be written in a 4 x 4-system for u^^, and 



p as: 

Pi — 
p 

Ptat _ 



VnVt + VnV„ + ^ = 0, (3 22) 

where we dropped B' from the system, since it is proportional to p' . Note that 7' vanishes. 
The system leads to the dispersion relation 

vt - c\l = 0, (3.23) 

which in differential form becomes: 

Apvlv', + vtp' - -fvlp't,, - 2^ptotVnv'„ - (2 - "f)BvlB' - (2 - -/)B^v„v'^ = 0. (3.24) 
Elimination of v'^, p' and B' gives 

^ - ^^3«2 + (^-2)c2^' • ^^-^^^ 
This expression allows us to then complete the exact solution as a function of 4>. 



4. Implementation and numerical details 

4.1. Details on the Newton- Raphson iteration 

We can generally note that ptot.pre < Ptot,post- This implies that the refraction has 3 
possible wave configurations: 2 shocks, a refiected rarefaction fan and a transmitted 
shock, or 2 expansion fans. Before starting the iteration on [[^]](p*), we determine the 
governing wave configuration. If [[^]](e) and [[^]](Ptot,5 — e) differ in sign, the solution 
has two rarefaction waves. If [[^]](ptot,5 + e) and [[^]](ptot,2 — e) differ in sign, the 
solution has a transmitted shock and a reflected rarefaction wave. In the other case, the 
solution contains two shocks in its configuration. If R is an expansion fan, we take the 
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Figure 4. The initial AMR grid at i = 0, for the example in section 5.1. 



guess 

min{ '——^ \ee {2,5}} +ptot,5 ^^^^ 

Po = — 2 ^ 

as a starting value of the iteration. This guess is the mean of the critical value ptot,crit , 
which satisfies 

«e,x - C^{Ptot,crit) = 0, (4.2) 

and P5, which is the minimal value for a transmitted shock. As we explain in section 5.3, 
'^2 x~ c^iPtot.crit) = is equivalent to ?;| — = and ^ — <?(ptot,crit) = is equivalent 
to V2 — <? = 0, and is thus a maximal value for the existence of a regular solution. If 
i? is a shock, we take (1 + e)ppost as a starting value for the iteration, where e is 10~^. 
We use a Newton-Raphson interation: p*^^ — p* — jr^J^, where f'{p*) is approximated 

numerically b 
where e = 10" 



numerically by ^"^"^ , where S = 10 ^. The iteration stops when ^' < e. 



4.2. Details on AMRVAC 
AMRVAC (|van der Hoist fc Keppe^ (|2007l ): iKeppens et~ai\ (|2003l )) is an AMR code, 



solving equations of the general form u* + V • F(u) = S(u,x, in any dimensionality. 
The applications cover multi-dimensional HD, MHD, up to special relativistic magneto- 
hydrodynamic computations. In regions of interests, the AMR code dynamically refines 
the grid. The initial grid of our simulation is shown in figure ID The refinement strategy 
is done by quantifying and comparing gradients. The AMR in AMRVAC is of a block- 
based nature, where every refined grid has 2^ children, and D is the dimensionality of 
the problem. Parallelisation is implemented, using MPI. In all the simulations we use 
5 refinement levels, starting with a resolution of 24 x 120 on the domain [0, 1] x [0,5], 
leading to an effective resolution of 384 x 1940. The shock is initially located at a; = 0.1, 
while the contact discontinuity is located at y — {x — \)tana. We used the fourth orde r 
Runge-Ku tta timestepping, together with a TVDLF-scheme fsee lToth fc Odstrcill (|l99d ): 



Ye3 (Il989l) 1 with Woodward-limiter on the primitive variables. The obtained numerical 



results were compared to and in agreement with simulations using other schemes, such 
as a Roe scheme and the TVD-Muscl scheme. The calculations were performed on 4 
processors. 

4.3. Following an interface numerically 

The AMRVAC implementation contains slight differences with the theoretical approach. 
Implementing the equations as we introduced them here would lead to excessive numerical 
diffusion on 7. Since 7 is a discrete variable we know 7(2;, y, t) exactly, if we are able to 
follow the contact discontinuity in time. Suppose thus that initially a surface, seperates 2 
regions with different values of 7. Define a function x : D x M+ — > IR. : (x, y, t) ^ x{^j 2/) ^)) 



a) 



2 4 6 



c) 



2 4 6 
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all shock solver - 



10 12 14 



no shock solver - 



10 12 14 



d) 
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2 4 6 



10 12 14 




Figure 5. [[:Jr"jj (P*) for the reference case from ISamtanevI 1120031 ): a) all shock solver; b) 
right shock solver; c) no shock solver; d) shock p* > pi. The all shock solver is selected. 



where D is the physical domain of {x, y). Writing x{^i u) = x{^j 0)j we ask x to vanish 
on the initial contact and to be a smooth function obeying 

• x{x,y) < 0, 

• 7 = 7r <^ x{x,y) > 0. 

We take in particular ±x to quantify the shortest distance from the point {x, y) to the 
initial contact, taking the sign into account. Now we only have to note that {xp)t = 
XPt+PXt ~ — • (/9v) — (pv- V)x = — V - (xpv). The implemented system is thus f l2.2p . 
but the last equation is replaced by {xpvx)^ + iXP'^y)y = 0. It is now straightforward to 
sh ow that we did not in troduce any new signal. In essence, this is the approach presented 
in lMulder 6^0!! (|l992l) . 



5. Results 

5.1. Fast-Slow example solution 

As a first hydrodynam ical examp l e, we set (a, 7;, 7^, ?y, M) — (|-, 0, |, |, 3, 2) , as 
originally presented in ISamtanevI (|2003l ). In figure [3 the first 3 plots show 
when assuming a prescribed wave configuration, for all 3 possible configurations. The 
last plot shows the actual function which consists of piecewise copies from the 

3 possible configurations in the previous plots. The initial guess is Pq = 4.111, the all 
shock solver is selected, and the iteration converges after 6 iterations with p* = 6.078. 
The full solution of the Riemann problem is shown in figure [H 
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Figure 6 . Solu tion to the fast-slowrefraction problem, for the reference case from ISamtanevI 
(l2003f ). Notice that p and — remain constant across the shocked contact. 




Figure 7. Solution to the slow-fast refraction problem from Ivan der Hoist fc Keppend (|2007l ). 
Notice that 5* remains constant across R. 
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Figure 8. Upper Left: p*{a). Note that for a < 0.61, there are no solutions for p*: the refraction 
is irregular; Upper Right: the wave pattern for regular refraction; Lower Left: For a = ^, the 
problem is 1-dimensional and there is no vorticity deposited on the interface. For decreasing a, 
the vorticity increases. Lower right: For regular refraction, {vy^^l > £5. 



5.2. Slow- Fast example 

In figure [7] we show the full solution of the HD Riemann problem, in which the reflected 
signal is an expansion f an, connected to the refract i on wit h parameters (a, , 7; , 7^ , ?y, M) 
(^,0,|,|,ji),10) from van der Hoist fc Keppend ( 2007 ). The refraction is slow-fast, and 



R is an expansion fan. Note that p and — remain constant across the CD, and the 
entropy S is an invariant across R. 

5.3. Tracing the critical angle for regular shock refraction 

Let us examine what the effect of the angle of incidence, a, is. Therefore we get back to the 
example from section 5.1, (/3~^,7/,7r.,r7, M) = (O, |,|,3,2) and let a vary: a G ] 0, . 
Note that a = f corresponds to a 1-dimensional Riemann problem. The results are shown 
in figure m Note that for regular refraction 5 > c|. We can understand this by noting 

that $± = '"■^■^"^■a^W^^-ge ^ ^ ^'e.x^e.yTCeV^e-Ce ^ ^ ^ whlch arc the eigenvalues 

of Gu"^ • Fu = (Fu^^ • Gu)-i. Note that we could have started our theory from the 
quasilinear form Uj, + (Gu^^-Fu)u3; =0 instead of equation p. lip . If we would have done 
so, we would have found eigenvalues ^, which would correspond to ^^^^^ . Moreover, both 
the eigenvalues, ^+ and have 4 singularities, namely £2 G {—Vx,2,Vx.2,—Vy,5,Vy^5} 
for ^_ and 65 € {—Vx,5,Vx,5,—Vy,2,Vy^2} for where thus C5 — vly ^ c\ = V2 and 
^2 = 2 = w| . It is now clear that it is one of the latter conditions that will be met 
for Oicrit- In the example, the transition to irregular refraction occurs at — Wy,5 — C5 and 
lim p* — '^^'■^^ (°'<:r-tt}-7i+i _ Q QY Figure [5] shows Schlieren plots for density 

from AMRVAC simulations for the reference case a = j, and the irregular case and 




Figure 9. Schlieren plots of the density for (/3 ^ ,'yi,'yr,ri, M) = (O, |, |, 3, 2) with varying a. 
Upper: a — a, regular reference case. Lower: a = 0.3: an irregular case. 



a = 0.3. In the regular case, all signals meet at the triple point, while for a < OLcrit — 0.61, 
the signals do not meet at one triple point, the triple point forms a more complex structure 
and becomes irregular. The CD, originated at the Much stem, reaches the triple point 
through an evanescent wave, which is visible by the contourlines. This pattern is called 
Mach Reflection- Refraction. Decreasi ng a even more, the ref lected wave transforms in 



a sequence of weak wavelets (see e.g. Nouragliev et al. ( 20051 )). This pattern, of which 



the case a — 0.3 is an example, is called Concave- Forwards irregular Refraction. These 
results are in agreement with our predictions. 

5.4. Abd-El-Fattah and Hendersons experiment 
In 1978, a shock tube experiment was performed bv lAbd-El-Fattah fc He nderson' (Il9786l). 
It became a typical test problem fo r simu lations (see e.g. lNouragliev et at. (2005i)) and re- 
fraction theory (see e.g. [Henderson ( 199ll )). The experiment concerns a slow-fast shock re- 
fraction at a C02/CH4 interface. The gas constants are 7CO2 = 1-288, 7c_ff4 = 1.303, /ico2 
44.01 and /xch^ = 16.04. Thus rj = = 0.3645. A very weak shock, M = 1.12 is 

refracted at the interface under various angles, von Neumann (Il943h theory predicts the 
critical angle acrit — 0.97 and the transition angle atrans — 1-01, where the reflected 
signal is irregular if a < Ucrit, a shock if acrit < a < atrans and an expansion fan if 
atrans < CX- This is in perfect agreement with the results of our solution strategy as illus- 
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1.38 




Figure 10. Exact solution for the Abd-El-Fattah experiment. Left: p*{a) confirms 
acrit = 0.97 and atrans — 1.01. Right: 0(a). 




0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.6 0.8 1 1.2 1.4 1.6 1.8 2 

Figure 11. Exact solution for (a, 7j, 7,-, M) = (-1,0,1,1,2) and a varying range of the 
density ratio r;. Left: for < 1 we have p* < ppoat = 4.5 and thus a reflected expansion fan, for 
> 1 we have p* > p^ost = 4.5 and thus a reflected shock. Right: for rj < 1: (jjT < ^ and for 
?7 > 1: (/>T > f . 

trated in figure [TOl There we show the pressure p* compared to the post shock pressure 
Ppost, as well as the angles (pR, 4>cd and 4>t for varying angle of incidence a. Irregular 
refraction means that not all signals meet at a single point. The transition at otcrit is 
one between a regular shock-shock pattern and an irregular Bound Precursor Refraction, 
where the transmitted signal is ahead of the shocked contact and moves along the contact 
at nearly the same velocity. This is also confirmed by AMRVAC simulations. If the angle 
of incidence, a, is decreased even further, the irregular pattern becomes a Free Precur- 
sor Refraction, where the transmitted signal moves faster than the shocked contact, and 
reflects itself, introducing a side- wave, connecting T to CD. When decreasing a even 
further, another transition to the Free Precursor von Neumann Refraction occurs. 

5.5. Connecting slow-fast to fast-slow refraction 

Another example of how to trace transitions by the use of our solver is done by changing 
the density ratio r] across the CD. Let us start from the example given in section 5.1 and 
let us vary the value of r]. 

Here we have (a, 7;, 7^, M) = (f ,0, |, |,2). The results are shown in figure [TTl 
Note that, since Ppost = 4.5, we have a reflected expansion fan for fast-slow refraction, 
and a reflected shock for slow-fast refraction. The transmitted signal plays a crucial 
role in the nature of the reflected signal: for fast-slow refraction 0t < f , but for slow- 
fast refraction, 0t > ■§ and the transmitted signal bends forwards. We ran our solver 
for varying values of M and a, and for all HD experiments with = jr, we came to 
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1,12 



rho 

14.67 
3.86 
3.04 
2.22 



Figure 12. Density plots for (q, 7;, 7^, M) = (-1,0,1,1,2). Left: A slow/fast refraction 
with Tj = 0.8. Note that > f and R is an expansion fan. Right: A fast/slow refraction with 
T) = 1.2. Note that 0t < f and i? is a shock. 




Figure 13. Left: Solution for the fast-slow problem: strong perpendicular magnetic fields de- 
crease the instability of the CD. Right: Solution for the slow-fast problem: strong perpendicular 
magnetic fields decrease the instability of the CD. 



the conclusion that a transition from fast-slow to slow-fast refraction, coincides with a 
transition from a reflected shock to a reflected expansion fan, with (/>t = §. This result 
agrees with AMRVAC simulations. In figure fT2l a density plot is shown for ?/ = 1.2 and 
r] = 0.8. 



5.6. Effect of a perpendicular magnetic field 

In general, the MHD equations result in the following jump conditions across a contact 
discontinuity 



(5.1) 



It follows, that if the component i?„ of the magnetic field, normal to the shock front is non- 
vanishing, a case we did not consider so far, the MHD equations do not allow for vortici ty 
deposition on a contact discontinuity and the RMI is suppressed ( Wheatlev et al. ( 2005f l). 
The remaining question is what the effect of a purely tangential magnetic field is, where 
the field is perpendicular to the shock front and thus acts to increase the total pressure 
and the according flux terms. 

Also note that it follows from equations 3.18 and 3.19 that is invariant across shocks 
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Figure 14. Density plots at t ^ 2.0 for (a, 7;, 7^, M) = (j,|,|,3, 2) with varying 
Upper: /J~^ = 0. The hydrodynamical Richtmyer-Meshkov instabihty causes the interface to 
roll up. Center: — i. Although the initial amount of vorticity deposited on the interface 
is smaller than in the HD case, the wall reflected signals pass the wall- vortex and interact with 
the CD, causing the RMI to appear. Lower: l3~^ — 1. The shock is very weak and the interface 
remains stable. 



and rarefaction fans. Therefore, can only jump across the shocked and unshocked 
contact discontinuity and B cannot change sign. 

Revisiting the example from section 5.1, we now let the magnetic field vary. Figure [T^ 
shows [[ft]](/?) across the CD. Also for rj — 0.8, making it a slow-fast problem, [[t't]](/3) is 
shown. First notice that no shocks are possible for /3 < 0.476, since u}+ would not satisfy 
0;+ > —M. Manipulating equation 13.91 we know that this is equivalent to 

P > Prmn = , (5.2) 
7;(M^ - 1) 

This relation is also equivalent to ci > M, which means that the shock is submagne- 
tosonic, compared to the pre-shock region. Figure [14] shows density plots from AMRVAC 
simulations at t = 2.0, for (a, 7;, 7^, ?7, M) = (|, |, |,3, 2) with varying f3~^. First note 
that the interface is instable for the HD case. Increasing (3~^ decreases the shock strength. 
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Figure 15. The reference problem from ISamtanevI l|2003l ) with varying l3. Left: The dependence 
of (pcD on /3. Note that ^ lim (f)CD = j = a, since this is the limit to infinitely weak shocks: 

lim At — Right: The vorticity deposition in the shocked contact scales as the Atwood 



number and lim 



= 1. 



For the interface remains stable, but for f3~^ = 1, the shock is very weak: the Atwood 
number At = 0.17, and the interface remains stable. 

Shown in figure [T51 is the vorticity across the CD. In the limit case of this minimal 
plasma-/? the interface is stable, both for fast-slow and slow-fast refraction. As expected, 
in the fast-slow case, the reflected signal is an expansion fan, while it is a shock in the 
fast-slow case. Also note that the signs of the vorticity differ, causing the interface to 
roll up clockwise in the slow-fast regime, and counterclockwise in the fast-slow regime. 
When decreasing the magnetic field, the vorticity on the interface increases in absolute 
value. This can be understood by noticing that the limit case of minimal plasma-/3 is also 
the limit case of very weak shocks. This can for example be understood by noting that 
^ lim (jjCD = ct (see figure [TSl ) . A convenient way to measure the strength of a shock 

is by use of its Atwood number 

At^P^^. (5.3) 

P2 + Pi 

Figure ITSl shows the jump across the shocked contact [[vt\], scaled to the shocks Atwood 
number. Note that in the limit case of very weak shocks the Atwood number equals the 
jump in tangential velocity across the CD, in dimensional notation: 



lim 



{[vt]] 



At 



(5.4) 



When keeping the Atwood number constant, the shocks sonic Mach number is given 



by 



M 



l + At / (2 - 27 - 7/3)Ai2 _^ + 2-i)At - -y^ ■ 



1-At 



(72/3)^^2 + (72/3 „ ji3)At - -f/3 



' (At + l)((7/3 + 27 - 2)At - (7/3 + 2)) 



y 7/3(1 -At) (7 At - 

Note that in the limit for weak shocks 



1) 



lim M 

At-^O 



hl3 + 2 
7/3 ' 



(5.5) 
(5.6) 

(5.7) 
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Figure 16. Left: Solution for the fast-slow problem: strong perpendicular magnetic fields de- 
crease the instability of the CD. Right: Solution for the slow-fast problem: strong perpendicular 
magnetic fields decrease the instability of the CD. 




Figure 17. AMRVAC plots of ^ for ylt = ^, with varying beta, upper: (3 = 16, lower: 

13 = 0.25. 

which is equivalent to l5.21 and in the limit for strong shocks, M — > oo. Figure [TBI shows 
the deposition of vorticity on the shocked contact, for a constant Atwood number. We 
conclude that under constant Atwood number, the effect of a perpendicular magnetic 
field is small: Stronger perpendicular magnetic field increase the deposition of vorticity 
on the shocked contact slightly. This is confirmed by AMRVAC sumulations (see figure 

[IZD. 



6. Conclusions 

We developed an exact Riemann solver-based solution strategy for shock refraction at 
an inclined density discontinuity. Our self-similar solutions agree with the early stages 
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of nonlinear AMRVAC simulations. We predict the critical angle acrit for regular refrac- 
tion, and the results fit with numerical and experimental results. Our solution strategy 
is complementary to von Neumann theory, and can be used to predict full solutions 
of refraction experiments, and we have shown various transitions possible through spe- 
cific parameter variations. For perpendicular fields, the stability of the contact decreases 
slightly with decreasing /3 under constant Atwood number. Wc will generalise our results 
for arbitrary uniform magnetic fields, where up to 7 signals arise. In this case we will 
search for non-evolutionary solutions, involving intermediate shocks, and for alternative 
evolutionary solutions, where the appearance of intermediate shocks can be avoided by 
including compound waves. We will investigate shock refraction involving initial slow, 
intermediate and fast shocks, and qualify the effect on the refraction. 
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